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1. Introduction 

Bessel function values I»(x) and J ,,{x) (x real), generated by BESLRI [l], 1 were compared 
with check values. For \x\< 64, these check values were calculated via the ascending series 

J " (x > { 2 ) £ o k\(n + k)V '" {X) \2) £ k\(n + k)V 

using multiprecision arithmetic [2]. For x ^ 64, the asymptotic formulas for large argument were 
used: 

/ / , e' ( Qx-1) (m-1)(/li-9) (m-1)( M -9)( M -25) 

V2^V 8 * 2 !(8a;) 2 3!(8*)'< 

where ju, = 4r-, and error is controlled by 7.16 of [3], and 

[Y 

J n {x) = \ — {P(n, x) cos X~ (?("> *) sin x} 

V 77% 



where x = x ~ I 9+7) 7r ' and 



p fn , , (m-1)(^-9) (m-1)(m-9)(m-25)( M -49) 
1 ' ; 2!(8x) 2 4!(8x) 4 

M-l (M-l)(M-9)(M-25) , 
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with errors bounded by the first neglected term [4, 9.2.9 and 9.2.10]. These were programmed in 
double precision, ror x ^ — 64, — x was used as the argument, and the sign of the function value 
was determined from the equations 

Jn(-X)=(-I)"j n (x), I II (-X)=(-1)»I„(X). 

The method used was that of bit comparison; the difference between the 60-bit mantissas of 
test and check values was expressed as a multiple, /n, say, of the last bit. Such a bit error cor- 
responds to a relative error between m • 2~~ m and m • 2~ 59 . 

This test is too strict near a zero of J n (x), where absolute error is a more realistic measure. 
Thus for n < \x\, both mantissas were right-shifted so that the bit error described above was a 
multiple of the 60th binary place. These cases are given above the dashed line in the computer 
printout; statistics for this test are compiled separately. 

A bit comparison test was also used to check the section of the code involving the two-term 
ascending series for small \x\. 

The error return feature was tested exactly as in [5]. 

2. The Bit Comparison Test 

To test the function I,,(x) , ten arguments were chosen in each interval 

Aj={x:2J~ 1 ^ \x\ <2->'}, ; = -13(l)10, 

the 60-bit mantissas of each x being produced by a pseudorandom number generator [5]. The orders 
used in all cases were n = 0(1)15. The only modification occurred in the set A\ {) , for which \x\ was 
limited to 700 to avoid computer overflow in the final answers. Summary statistics for the 60-bit 
test on I n are: 

373 results correct to 60 bits 

583 results in error by 1 multiple of the last bit 

383 results in error by 2 multiples of the last bit 

288 results in error by 3 multiples of the last bit 

232 results in error by 4 multiples of the last bit 

196 results in error by 5 multiples of the last bit 

176 results in error by 6 multiples of the last bit 

130 results in error by 7 multiples of the last bit 

1479 results in error by more than 7 multiples of the last bit 



3840 (total) 

The largest bit error was (36)s, implying that at least 55 bits are correct in each case. 

For the function J n (x), the set of intervals Aj was expanded toj=— 13(1)16. Again, ten pseudo- 
random values were used in each interval, and the orders used were 0(1)15. Summaries for the 
tests are: 

Test on 60 Significant Bits \x\ < 64 

403 results correct to 60 significant bits 

643 results in error by 1 multiple of the last bit 

453 results in error by 2 multiples of the last bit 

331 results in error by 3 multiples of the last bit 
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232 results in error by 4 multiples of the last bit 

171 results in error by 5 multiples of the last bit 

109 results in error by 6 multiples of the last bit 

68 results in error by 7 multiples of the last bit 

82 results in error by more than 7 multiples of the last bit 



2492 (total) 

The largest bit error was (15)«, showing that at least 56 bits are correet in each ease. 

Test on 60 Binary Places \x\ < 64 

161 results correct to 60 places 

200 results in error by 1 multiple of the last place 

100 results in error by 2 multiples of the last place 

85 results in error by 3 multiples of the last place 

56 results in error by 4 multiples of the last place 

46 results in error by 5 multiples of the last place 

45 results in error by 6 multiples of the last place 

14 results in error by 7 multiples of the last place 
1 result in error by more than 7 multiples of the last place. 

708 (total) 

The largest bit error was (10)*, showing that 56 binary places are correct in each case, and 57 
places in all but one case. 

In the test on 60 binary places for \x\ ^ 64, the largest error was (422)g. This occurred in ./:,(.%) 
and J 7 (x) for x = 59766.6 .... Bearing in mind that about 60,000 back-recursion steps were 
executed, that each step included two multiplications and one addition, and also that the computer 
uses truncation in floating-point arithmetic rather than roundoff, we see that the error of (422) s 
= 274 is quite reasonable. The other errors can be explained similarly. 

For the test of the two-term ascending series, arguments of 2~ Ui and 2~ 100 were used; the former 
is just barely small enough for the truncation error to be neglected. Both I n (x) and J „{x) were 
tested with n= 0(1) 15, and the largest error was 6 multiples of the last bit. 

3. Testing of Error Return 

The method used is the same as in [5 J. Arguments used were x = 2 k where k = —13(1)9 for 
/'s, and /c = — 13(1)13 for J's. In the test on Bessel function values of order NCALC-1, the largest 
error was 6 multiples of the last bit; most were zero. 

In the test on order [x], the largest error was (46)8, which occurred in 7hi92(8192). This test 
value involved about 1600 more back-recursion steps than were required for the check value, so 
this error is quite reasonable. The other errors were (35)s for jc = 4096, zero for x = 2048, (12)« 
for x= 1024; for x ^ 512, the errors did not exceed 5. 

4. Summary and Conclusions 

For n = 0(l)15 and the values of x tested, BESLRI yielded J „(x) and I n (x) with either a rela- 
tive error or an absolute error bounded by 10~ 14 . The latter assessment applied to J „{x) when 
n< \x\. For \x\ < 64, the bound was 10~ 16 . The error return and ascending series features were 
tested and found to be executed accurately in all cases. 
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For other values of n and x, similar accuracy can be expected, since the tests we have used 

|*| 
include large and small values of \x\ and , the former as large as 2 16 . 

Calculations were performed on a UNIVAC 1108 under EXEC 2. All test values were cal- 
culated with double precision mantissa of 60 bits, corresponding to just over 18 figures. Subroutine 
executions took approximately 0.2 N milliseconds, where 

A=max (|*|,NB) + 10. 

5. Appendix: Algorithm BESLRI 

a ELT BESLCI, 1,7301 12, 5174a 
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SUBROUTINE BESLCH X,Y,NB, IZE,BR, BI.NCALC ) 
THIS ROUTINE CALCULATES BESSEL FUNCTIONS I AND 
COMPLEX ARGUMENT AND INTEGER ORDER. 



EXPLANATION OF VARIABLES IN TEE CALLING SEQUENCE 

DOUBLE PRECISION REAL PART OF THE COMPLEX ARGUMENT 
FOR WHICH I*S OR J*S ARE TO BE CALCULATED. IF I*S 
ARE TO BE CALCULATED, ABS{ X ) MUST NOT EXCEED EXPARG 
(WHICH SEE BELOW ). 

IMAGINARY PART OF THE ARGUMENT. IF J»S ARE TO BE 
CALCULATED, ABS(Y) MUST NOT EXCEED EXPARG. 
INTEGER TYPE. 1 ♦ HIGHEST ORDER TO BE CALCULATED. 
IT MUST BE POSITIVE. 

INTEGER TYPE. ZERO IF J*S ARE TO BE CALCULATED, 1 
IF I*S ARE TO BE CALCULATED. 

DOUBLE PRECISION VECTOR OF LENGTH NB, NEED NOT BE 
INITIALIZED BY USER. IF THE ROUTINE TERMINATES 
NORMALLY, (NCALC-NB), IT RETURNS THE REAL PART OF 
J( OR I)-SUB-ZERO THROUGH J(0R I )-SUB-NB-MINUS -ONE 
OF Z IN THIS VECTOR. 
IMAGINARY ANALOG OF BR. 
C INTEGER TYPE, NEED NOT BE INITIALIZED BY USER. 

BEFORE USING THE RESULTS, THE USER SHOULD CHECK THAT 
NCALC-NB, I.E. ALL ORDERS HAVE BEEN CALCULATED TO 
THE DESIRED ACCURACY. SEE ERROR RETURNS BELOW. 



C EXPLANATION OF MACHINE- DEPENDENT CONSTANTS 

C 

C NSIG DECIMAL SIGNIFICANCE DESIRED. SHOULD BE SET TO 

C IFIX( ALOG10(2 )#NBIT*-l ), WHERE NBIT IS THE NUMBER OF 

C BITS IN THE MANTISSA OF A DOUBLE PRECISION VARIABLE, 

C SETTING NSIG HIGHER WILL INCREASE CPU TIME WITHOUT 

C INCREASING ACCURACY, WHILE SETTING NSIG LOWER WILL 

C DECREASE ACCURACY. IF ONLY SINGLE-PRECISION 

C ACCURACY IS DESIRED, REPLACE NBIT BY THE NUMBER OF 

C • BITS IN THE MANTISSA OF A SINGLE-PRECISION VARIABLE. 

C THE RELATIVE TRUNCATION ERROR IS LIMITED TO T-.5*10 

C *»-NSIG FOR ORDER GREATER THAN ABS(Z), AND FOR ORDER 

C LESS THAN ABS(Z) (GENERAL TEST), THE RELATIVE ERROR 

C IS LIMITED TO T FOR FUNCTION VALUES OF MAGNITUDE AT 

C LEAST 1, AND THE ABSOLUTE ERROR IS LIMITED TO T FOR 

C SMALLER VALUES. 

C NTEN LARGEST INTEGER K SUCH THAT 10*#K IS MACHINE- 

C REPRESENTABLE IN DOUBLE PRECISION. 

C LARGEZ UPPER LIMIT ON TEE MAGNITUDE OF Z. BEAR IN MIND 

C THAT IF ABS(Z)-N, THEN AT LEAST N ITERATIONS OF THE 

C BACKWARD RECURSION WILL BE EXECUTED. 

C EXPARG LARGEST DOUBLE PRECISION ARGUMENT THAT THE LIBRARY 

C DEXP ROUTINE CAN HANDLE. 
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000055 C ERROR RETURNS 00550 

000056 C 00560 

000057 C LET G DENOTE EITHER I OR J. 00570 
00005e C IN CASE OF AN ERROR, NCALC.NE.NB, AND NOT ALL G»S 00580 

000059 C ARE CALCULATED TO THE DESIRED ACCURACY. 00590 

000060 C IF NCALC. I.T.O, AN ARGUMENT IS OUT OF RANGE. NB.LE.O 00600 

000061 C OR IZE IS NEITHER NOR 1 OR IZE-0 AND ABS( Y ) .GT .E XPARG , 00610 

000062 C OR IZEM AND ABS( X ). GT. EXP ARC . IN THIS CASE, TEE VECTORS 00620 

000063 C BR AND BI ARE NOT CALCULATED, AND NCALC IS SET TO 00630 
00C064 C MIN0(NB,0)-1 SO NCALC.NE.NB. 00640 

000065 C • NB. GT. NCALC. GT. WILL OCCUR IF NB.GT.MAG7. AND A3S(G- 00650 

000066 C SU3-NB-OF-Z/G-SUB-MAGX-OF-Z ).LT. 10.**< NTEN/2 )» I. E. NB 00660 
CO0067 C IS MUCH GREATER THAN MAGZ. IN THIS CASE, BR(N) AND BI(N) 00670 

000068 C ARE CALCULATED TO THE DESIRED ACCURACY FOR N . LE . NCALC, 00680 

000069 C BUT FOR NCALC. LT. N. LE. NB, PRECISION IS LOST. IF N.GT. 00690 

000070 C NCALC AND ABS( G( NCALC- 1 )/G( N - 1 ) ) . EQ. I 0« * - K, THEN THE LAST 00700 

000071 C K SIGNIFICANT FIGURES OF G(N-l) ( «BR( N ) ♦ I *B I ( N ) ) ARE 00710 

000072 C ERRONEOUS. IF THE USER WISHES TO CALCULATE G(N-l) TO 00720 
00 0073 C HIGHER ACCURACY, HE SHOULD USE AN ASYMPTOTIC FORMULA FOR 730 
00007* C LARGE ORDER. 00740 

000075 C 00750 

000076 DOUBLE PRECISION 00760 

000077 1 X,Y,BR,BI,PR,PI, PLASTR,PLAST I, POLDR, POLD I .PSAVER, OOV/0 

000078 2 PSAVEI,EXPARG,TEST,TOVER,TEtfPAR,TEMPAI,TEMPBR,TEMPBI, 00780 

000079 3 TEMPCR, TEMPCI , SIGN, SUMR, SUMI , ZINVR, ZINVI 00790 
000030 DIMENSION BR< NB ), BI( NB ) 00800 
00008L DATA NS IG, NTEN, LARCEZ , EXPARG/ 19, 307, 1 000, 7 ,D2/ 00810 
000082 TEMPAR-DSQRTC X#X»Y»Y ) 00820 
OC0083 MAGZ-IFIX( SNGL( TEMPAR )) 00830 

000084 IF( NB.GT.O. AND. MAGZ. LE. LARGEZ. AND. (( IZE. EQ.O. AND. 00840 

000085 1 i>ABS< Y ).LE.EXPARG ).OR. ( IZE.EQ.l ,AND.DAES( X ) . LE . 00850 

000086 2 EXPARG))) GO TO 1 00860 

000087 C ERROR RETURN -- Z, NB, OR IZE IS OUT OF RANGE 00870 

000088 NCALC"MINO( NB, )- 1 00880 

000089 RETURN 00890 

000090 1 SIGN-DBLEC FLOATC 1-2»IZE )) 00900 

000091 NCALC-NB 00910 

000092 C USE 2 -TERM ASCENDING SERIES FOR SMALL Z 00920 

000093 IF( TEMPAR**4. LT. . 1D0#«NSIG ) GO TO 50 00930 
00 0094 C INITIALIZE THE CALCULATION (iF THE P»S 00940 

000095 NBMZ-NB-MAGZ 00950 

000096 N-UAGZM 00960 

000097 IF( DABS( X ). LT. DABS( Y ) ) GO TO 2 00970 
OC0O93 ZINVR-1 .D0/( X*Y»Y/X ) 00980 

000099 ZINVI--Y*ZINVR/X 00990 

0001 00 GO TO 3 1 000 

000101 2 ZINVI --1 .D0/( Y*X«X/Y) 01010 

000102 ZINVR--X»ZINVI/Y 1020 

000103 3 PLASTR-1.D0 1030 
000104- PLASYI-O.DO 01040 

000105 PR-SIGN*DBLE( FLOATC 2»N) )*ZINVR 01050 

000106 PI-SICN*D3LE< FLOATC 2*N ) )»ZINVI 01060 
00010/ • TEST-2.D0»1 .D1*»NSIG 01070 

000108 M-0 01080 

000109 IF( NBMZ.LT.3 ) GO TO 6 01090 

000110 C CALCULATE P*S UNTIL N-NB-1. CHECK FOR POSSIBLE OVERFLOW. O1100 

000111 TOVER-1 .Dl*»( NTEN-NSIG) OHIO 

000112 NSTART-MAGZ*2 0X120 
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NEND-NB-1 

DO 5 N«NSTART,NEND 

POLDR -PLASTR 

POLDI-PLASTI 

PLASTR-PR 

PLASTI-PI 

PR-SIGN#( DBLE<FL0AT( 2#N ) )#( PLASTR«ZINVR-PLASTI#ZINVI 
1 )-P0LDR) 

PI-SIGN»< DBLE( FLOAT( 2*N ) )»( PLASTI*ZINVR*PLASTR*ZINVI 
1 )-POLDI) 

IF( ( PR/T0VER)*»2MPI/T0VER)#*2-1 .DO ) 5,5,7 

5 CONTINUE 
N-NEND 

C CALCULATE SPECIAL SIGNIFICANCE TEST FOR NBMZ.GT.2. 
TEMPBI-DMAXH DABS(PR),DABS( PI ) ) 

TEMPBI-TEMPBI*DSQRT(2.D0*1.D1*#NSIG#DSQRT( ( ( PR/TEMPBI ) 
1#*2 ♦( PI/TEMPBI )**2 )*( (PLASTR/TEMPBI )**2 ♦( PLASTI / 
2TEMPBI )**2))) 
TEST-DWAXK TEST.TEMPBI ) 
C CALCULATE P» S UNTIL SIGNIFICANCE TEST IS PASSED. 

6 N-N*l 
POLDR-PLASTR 
POLDI-PLASTI 
PLASTR-PR 
PLASTI-PI 

PR -SIGN*( DBLE(FL0AT( 2*N ) )*( PLASTR»ZIN VR-PLASTI*ZINVI ) 
1 -PeLDR) 

PI-SIGN*( DBLE(FLOAT( 2*N ) )#( PLASTI«ZIN VR*PLASTR*ZINVI ) 
1 -POLDI) 

IF( ( PR/TEST )**2*( PI/TEST)*#2.LT. 1. DO) GO T6 6 

IF( M.EO. 1 ) GO TO 12 
C CALCULATE STRICT VARIANT OF SIGNIFICANCE TEST, AND 
C CALCULATE T*S UNTIL THIS TEST IS PASSED. 

M-l 

TEMPUI-DMAXK DABS( PR),DABS( PI ) ) 

TEUPBR-DSORT( ( (PR/TEMPBI )»»2*( P I/TEMPB I )»»2)/ 
1 ( ( PLASTR/TEMPBI )#*2*( PLASTI/TEMPBI )»*2 ) ) 

TEMPBI«DBLE(FLGAT(N»1 ) )/TEMPAR 

IFC TEMPBS* 1 „D0/TEMPBR.GT.2.D0*TEMPBI ) TEMPBR-TEMPBI 
1 ♦D£CRT( TEMPBI«-*2-l .DO ) 

TEST-TEST/DSQRT(TEMPBR-1.D0/TEMPBR) 

IF( (PR/TEST>»*2*( PI/TEST)*»2-1.D0 ) 6, 12; 12 

7 NSTART-N*1 

C TO AVOID OVERFLOW, NORMALIZE P»S BY DIVIDING BY TOVER. 
C CALCULATE P* S UNTIL UNN0RUALIZED P WOULD OVERFLOW. 

PR -P2/ TOVER 

PI-PI/T0VER 

PLASTR - PLASTR /TOVER 

PL ASTI - PLASTI /TOVER 

P SAVER -PR 

PSAVEI-PI 

TEMFCR- PLASTR 

TEMPCI- PLASTI 

TEST-1 . D1**(2*NSIG ) 

8 N-N*l 

POLDR -PLASTR 
POLDI-PLASTI 
PLASTR-PR 
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0001 7 1 
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17 3 
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000201 
000202 
000203 
000204 
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000207 
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000210 
00021 t 
000212 
000213 
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00021 5 
000216 
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PLASTI -PI 

PR-SIGN«( DBLH( FLOAT( 2 *N ) )*( PLASTR*ZI N V tf- PLASTI* ZINVI ) 
1 -I'Ot.nR) 

PI-S1GN*< DBLE( FLOAT( 2»N ) )*< PLASTI «ZIN V K *PLASTR*ZIN VI ) 
1 -POLDI ) 
IF( PR**2*PI**2.LE. TEST) GO TO 8 
C CALCULATE BACKWARD TEST, AND FIND NCALC, THE HIGHEST N 
C SUCH THAT THE TEST IS PASSED. 

TEMPBR-D£QRT( ( PLASTR**2*PLASTI**2 )/( PO LDR*«2 ♦P0LDI**2 ) 
1 ) 
TEMPBI -DBLE( FL0AT( N ) )/TEMPAR 

IPC TEMPBRM . D0/TEMPBR.GT. 2.D0*TEMPBI ) TEMPER -TEMPB I • 
I DSQRT( TEMPBI**2- 1 . DO ) 
TEST- . 5D0*( 1 . DO- 1 . DO/TEMP BR*»2 )/l .D1**NSIG 
TEST-( ( PLASTR»*2*PLASTI**2 )*TEST)*( ( P0LDR**2 »POLDI*»2 ) 
1 *TEST) 
PR-PLASTR»TOVER 
PI-PLASTI*TOVER 
N-N-l 

NEND-MINO( NB,N) 
DO 9 NCALC- NSTART. NEND 
POLDR-TEMPCR 
P0LDI-TEMPCI 
TEMPCR-PSAVER 
TEMPCI-PSAVEI 

PSAVER«SIGN*(DBLE(FLOAT( 2*N ) )*( TEMPC R*ZINVR-TEMPCI» 
1 ZINVI )-POLDR ) 

PSAVEI-SIGN*( DBLE( FL0AT( 2*N ) )*( TEMPC I *ZIN VR*TEMPCR* 
1 ZINVI )-POLDI ) 

IF( ( PSAVER**2*PSAVEI**2 )»( TEMPC R**2 *TEMPCI»*2 )-TEST) 
1 9,9, 10 
9 CONTINUE 

NCALC-NEND* 1 
10 NCALC-NCALC-1 
C THE COEFFICIENT OF B( N ) IN THE NORMALIZATION SUM IS 
C M*SQRT( -1 )**IMAG, WHERE M--2.0, OR 2, AND I MAG IS OR 1. 
C CALCULATE RECURSION RULES FOR M AND IMAG, AND INITIALIZE 
C THEM. 

12 N-N*l 

TEMPBR-DBLE(FL«AT( IZE ) )*X»DBLE( FL0ATC 1 -IZE ) )»Y 

IPOS-0 

IF(TEMPBR) 13,14,13 

13 IPOS-IFIX< SNGL( 1 . 1 DO* TEMP BR/ DA BS( TEMPB R) >) 

14 MRhCUR«4»( ( 2* IZE* IPOS )/2)-3-2*( IZE* IPOS) 
K-2*IP0S*2* IZE*riPOS««2- IZE 

L'N-4»( N/4 ) 

MLAST-2*«*( ( K»L>/4 )-4»( < K*L)/2 ) 

IF( IPOS.EQ. 0. AND.( L.EQ. 1. OR. L.EQ. 3}) MLAST-0 

L-L*3-4*( ( L*3 )/4 ) 

M -2*8*< < K*L)/4 )-4*( ( K»L)/2 ) 

IF( IPOS.EQ. 0. AND.( L.EQ. 1.0R.L.EO.3)) M-0 

IMRECR-( 1-IZE )»IPOS**2 

IMAG- IMRECR*( L-2*( L/2 ) ) 
C INITIALIZE THE BACKWARD RECURSION AND THE NORMALIZATION 
C SUM. 

TEMPBR-O.DO 

TEMPBI-O.DO 

IF( DABS( PI ).GT.DABS( PR) ) GO TO 15 
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000274 
000275 
000276 
000277 
000278 
000279 
000280 
000281 
000282 
0002e3 
000284 
000285 
000286 



TEMPAR-1 .D0/( PR*PI»( PI/PR) ) 
TEMPAI--( PI*TEMPAR)/PR 
Gfl T6 16 

15 TEMPAI--l.D0/( PI*PR*( PR/PI ) ) 
TEMPAR-.( PR*TEMPAI )/PI 

16 IF( IMAG.NE.O ) GO TO 17 
SUMR«DBLE(FLOAT(M) )»TEMPAR 
SUMI-DBLE(FLOAT(M) >»TEMPAI 
GO TO 18 

17 SUMR--DBLE( FL0AT( M ) )» TEMP A I 
SUMI-DBLE(FLOAT(M) )*TEMPAR 

18 NEND-N-NB 
IF(NEND) 26,22,19 

C RECUR BACKWARD VIA DIFFERENCE EQUATION CALCULATING (BUT 
C NOT STORING) BR( N ) AND BI(N) UNTIL N-NB. 

19 DO 21 L-1,NEND 

N-N-l 

TEMPCR-TEMPBR 

TEMPCI-TEMPBI 

TEMPBR-TEMPAR 

TEMPBI-TEMPAI 

PR-DBLE( FLflATC 2*N ) )*ZINVR 

PI-DBLE( FL0AT(2#N))»ZINVI 

TEMPAR-PR*TEMPBR-PI»TEMPBI-SIGN»TEMPCR 

TEMPAI-PR*TEMPBI*PI*TEMPBR-SIGN»TEMPCI 

I MAG -( 1 - I MAG )« I MRECR 

K -ML AST 

MLAST-M 

M-IOMRECUR 

IF( IMAG.NE. ) GO TO 20 

SUMR-SUMR*DBLE(FLOAT( M))»TEMPAR 

SUMI-SUMI*DBLE(FLOAT( M))« TEMPAI 

GO TO 21 

20 SUMR-SUMR-DBLEC FLOAT( M) )*TEMPAI 
SUMI-SUMI*DBLE( FLOAT( M) )«TEMPAR 

21 CONTINUE 

C STORE BR(NB), BI( NB ) 

22 BR(N)-TEMPAR 
BI( N)-TEMPAI 

IF( N.GT.l ) GO TO 23 
C NB-1. SINCE 2*TEMPAR AND 2*TEMPAI WERE ADDED TO SUMR AND 
C SUM I RESPECTIVELY, WE MUST SUBTRACT TEMPAR AND TEMPAI 

SUMR-SUMR-TEMPAR 

SUM I- SUM I -TEMPAI 

G6 TO 35 
C CALCULATE AND STORE BR( NB- 1 ), BI ( NB- 1 ) 

23 N-N-l 

PP-DBLE( FL0AT(2*N ) )»Z INVR 

PI-DBLE(FL0AT(2*N) )»ZINVI 

BR( N)-PR*TEMPAR-PI*TEMPAI-SIGN#TEMPBR 

BI( N)-PR*TEMPAI*PI*TEMPAR-SIGN*TEMPBI 

IF( N.EO. 1 ) G6 TO 34 

IMAG-C 1-IMAG)*IMRECR 

K-MLAST 

MLAST-M 

M-K*MRECUR 

IF( IMAG.NE.O ) GO TO 24 

SUMR-SUMR«DBLE( FLOAT( M) )»BR(N ) 



02290 
02300 
02310 
02320 
02330 
02340 
02350 
02360 
02370 
02380 
02390 
02400 
02410 
02420 
02430 
02440 
02450 
02460 
02470 
02480 
02490 
02500 
02510 
02520 
02530 
02540 
02550 
02560 
02570 
02580 
02590 
02600 
02610 
02620 
02630 
02640 
02650 
02660 
02670 
02680 
02t>90 
02700 
02710 
02720 
02730 
02740 
02750 
02760 
02770 
02780 
02790 
02800 
02810 
02820 
02830 
02840 
02850 
02860 
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0OOP87 
000288 
000289 
000290 
0O0291 
000292 
000293 
000294- 
000295 
000296 
000297 
000298 
000299 
0300 
000301 
000302 
000303 
000304 
000305 
000306 
000307 
000308 
0003C9 
000310 

00031 1 

00031 2 
000313 
0314 
0C0315 
000316 
000317 
0C0318 
000319 
000320 
000321 
000322 
000323 
000324 
000325 
000326 
000327 
000323 
0C0329 
000330 
000331 
000332 
000333 
000334 
000335 
000336 
000337 
000338 
000339 
000340 
000341 
000342 
000343 
000344 



000345 
000346 
000347 
000348 
000349 
O0C35O 
000351 
000352 
000353 
000354 
000355 
0C0356 
000357 
000358 
000359 
000360 
00036 1 
000362 
000363 
000364 
000365 
000366 
000367 
000368 
000369 



SUMI-SUMI*DBLEC FLOATC M) )«BI( N) 
G0 TO 30 
24 SUMR-SUMR-DBLEC FLOATC M) )*BI( N) 
SUMI-SUMI*DBLEC FLOATC M) )»BRC N > 
G6 TO 30 
C N.LT.NB, S6 STORE BR( N ), BI(N), AND SET HIGHER ORDERS ZERO 

26 BR( N)-TEMPAR 
BI( N)-TEMPAI 
NEND--NEND 

DO 27 L-l.NEND 
BR( N»L)-0. DO 

27 BI( N*L)-0.D0 

30 NEND-N-2 

IF( NEND.EQ.O ) GO TO 33 
C CALCULATE VIA DIFFERENCE EQUATION AND STORE BR(N).BKN), 
C UNTIL N-2 

D0 32 L-l.NEND 

N-N-l 

PR-DBLEC FLOATC 2*N) )*ZINVR 

PI-DBLE( FLOAT( 2*N) >»ZINVI 

BR( N )-PR*BR(N*l )-PI»BI( N*l )-SIGN*BR( N»2 ) 

BI( N )-PR»BIC N*l )*PI»BRC N*l )-SIGN*BI( N*2 ) 

IMAG"( 1-IMAG >*IMRECR 

K-MLAST 

MLAST-M 

M-K»MRFCUR 

IF( IMAG.NE. ) GO TO 31 

SUMR-SUMR*DBLEC FLOAT( M) )#BRC N ) 

SUMI-SUMI*DBLEC FL0AT( M> )«BI( N ) 

GO TO 32 

31 SUMR-SUMR-DBLEC FLOATC M) )»BI( N ) 
SUMI-SUMI»DBLE( FLOAT( M) )*BR( N ) 

32 CONTINUE 

C CALCULATE AND STORE BR( 1 ), BI(1) 

33 BR( 1)-2.D0»(HR(2)*ZINVR-BI(2)*ZINVI )-SIGN*BRC 3) 
BI( 1 )-2.D0*( IJR( 2 )»ZINVI*BI( 2)»ZINVR >-SIGN*BIC 3) 

34 SUMR-SUMR»BR( 1 ) 
SUMI-SUMI»8IC 1 ) 

C CALCULATE NORMALIZATION FACTOR, TEMPAR ♦I«TEMPAI 
3 5 IF(IZE.EO.l) GO TO 36 

TEMPCR-DBLEC FLOATC IPOS) )»Y 
TEMPCI-DBLE( FLOATC - IPOS ) )»X 
GO TO 37 . 

36 tempcr-dblec floatc ipos) )»x 
tempci-dblec float( ipos) )*y 

37 tempcr-dexpc tempcr ) . , 
tempbr-dcosc tempci ) 

tempbi-dsinc tempci ) 

ifc dabsc sumr ).lt. dabsc sumi ) ) go to 38 

tem pc i - su m i / su mr 

tempcr-c tempcr/sumr )/( 1 . do *tempc i*tempci ) 

tempar- tempcr»ctempbr*tempbi*tempci ) 

tempai-7empcr*( tempbi -temp br#tempc i ) 

GO TO 39 

38 TEMPCI-SUMR/SUMI 

TEMPCR-C TEMFCR/SUMI )/( 1 . DO ♦ TEMPCI »TEMPCI ) 
TEMPAR- TEMPCR*C TEMPBR*TEMPC I*TEMPBI J 
TEM PA I - TEMPCR* ( TEMPBI * TEMPC I - TEMPBR ) 



C NORMALIZE 

39 DO 40 N-l.NB 

TEMPBR-BRCN )»TEMPAK-B I( N )*TEMPAI 
BI( N )-BR( N)*TEMPAI*BI(N )*TEMPAR 

40 BR( N ) -TEMPBR 
RETURN 

C TWO-TERM ASCENDING SERIES FOR SMALL Z 

50 TEMPAR- I. DO 
TEM PA I -0. DO 

TEMPCR -.2 5 DOM X*X-Y#Y ) 

TEMPCI-.5D0*X*Y 

BR( 1 )-l .DO-SIGN*TEMPCR 

BI( 1 )--SIGN*TEMPCI 

IF(NB.EQ.l) GO TO 52 

DO 51 N-2.NB 

TEMPBR-C TEMPAR»X-TEMPAI*Y )/DBLE( FLOAT( 2*N-2 ) ) 

TEMPAI-( TEMPAR»Y*TEMPAI»X )/DBLEC FLOATC 2#N-2 )) 

TEMPAR- TEMPBR 

TEMPBR-DBLEC FLOATC N ) ) 

BRC N)-TEMPAS«C 1 . DO- SIGN*TEMPCR/TEKPB R ) 
1 ♦TEMPAKTEMPCI/TEMPBR 

51 BIC N )-TEMPAI*C 1 . DO- SIGN»TEMPCR/TEMPBR ) 
1 -TEMPAR*TEMPCI/TEMPBR 

52 RETURN 
END 



END CUR LCC 1102-0038 L8 



02870 
•02880 
02890 
02900 
02910 
02920 
02930 
02940 
02950 
02960 
02970 
02980 
02990 
03000 
03010 
03020 
03030 
03040 
03050 
03060 
03070 
03000 
03090 
03100 
031 10 
03120 
03130 
03140 
03150 
031&0 
03170 
03180 
03190 
03200 
03210 
03220 
03230 
03240 
03250 
03260 
03270 
03280 
03290 
03300 
03310 
03320 
03330 
03340 
03350 
03360 
03370 
03380 
03390 
03400 
03410 
03420 
03430 
,03440 



03450 
03460 
03470 
03480 
03490 
03500 
03510 
03520 
03530 
03540 
03550 
03560 
03570 
03580 
03590 
03600 
03610 
03620 
03630 
03640 
03650 
03660 
03670 
03680 
03690 
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